Translation of collision geometry fluctuations into momentum anisotropics in 

relativistic heavy-ion collisions 



O 

(N 

00 



Guang-You Qin, Hannah Petersen, Steffen A. Bass and Bcrndt Miiller 
Department of Physics, Duke University, Durham, NC 27708, USA 
(Dated: September 21, 2010) 

We develop a systematic framework for the study of the initial collision geometry fluctuations in 
relativistic heavy-ion collisions and investigate how they evolve through different stages of the fireball 
history and translate into final particle momentum anisotropics. We find in our event-by-event 
analysis that only the few lowest momentum anisotropy parameters survive after the hydrodynamical 
evolution of the system. The geometry of the produced medium is found to be affected by the pre- 
cquilibrium evolution of the medium and the thermal smearing of the discretized event-by-event 
initial conditions, both of which tend to smear out the spatial anisotropies. We find such effects to 
be more prominent for higher moments than for lower moments. The correlations between odd and 
even spatial anisotropy parameters during the pre-equilibrium expansion are quantitatively studied 
and found to be small. Our study provides a theoretical foundation for the understanding of initial 
state fluctuations and the collective expansion dynamics in relativistic heavy-ion collisions. 
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I. INTRODUCTION 

Experiments at the Relativistic Heavy Ion Collider 
(RHIC) have discovered that the strongly interacting 
matter produced in these highly energetic collisions ex- 
hibits strong collective flow, which can be well described 
by relativistic hydrodynamics ^ n noncentral colli- 

sions, the collective flow is azimuthally asymmetric in the 
plane transverse to the beam axis. This has been under- 
stood as the consequence of the initial spatial asymmetry 
of the medium produced by the two colliding nuclei which 
translates into a momentum anisotropy of the emitted 
particles due to the hydrodynamic expansion of the mat- 
ter. The magnitude of this flow anisotropy is quantified 
by the Fourier expansion coefficients v n of the azimuthal 
angular distribution of the emitted particles in the trans- 
verse plane @. 

The elliptic flow V2 signal has been extensively stud- 
ied in Au+Au collisions at RHIC as a function of various 
quantities (Toj . Hydrodynamic simulations have shown 
that elliptic flow V2 is sensitive to various transport prop- 
erties of the expanding hot medium, especially the spe- 
cific shear viscosity 77, the presence of which tends to re- 
duce the amount of the elliptic flow that can be built up 
in an ideal hydrodynamical fluid fllT - flsT ]. Considerable 
effort has been devoted to the quantitative extraction 
of the shear viscosity by comparing the measured elliptic 
flow V2 with viscous relativistic hydrodynamic simulation 
of the fireball evolution and other Boltzmann transport 
models that involve the violation of ideal hydrodynamic 
behavior flrl These comparisons have yielded an 

upper limit for the shear viscosity to entropy density 
s ratio: r\j s < 0.5 Ojjl, the same order of magni- 
tude as the conjectured KSS bound rj/s = l/(47r) [201 ] . 
which were obtained using anti-dc-Sittcr/conformal field 



theory (AdS/CFT) correspondence for certain quantum 
field theories similar to QCD. 

Current efforts in the extraction of the shear viscos- 
ity from precise v-i measurements arc subjected to vari- 
ous uncertainties in the hydrodynamic simulations, i.e., 
equations of state, large shear viscosity in late hadronic 
stage [H]], bulk viscosity [13], and the treatment of 
the freeze-out conditions. Among the largest uncertain- 
ties is the initial geometry employed in the hydrody- 
namical simulations, i.e^jthc initial fireball eccentricity 
£2 = (y 2 — x 2 ) I '(y 2 + x 2 ) 0,|2j|. In ideal hydrodynamics, 
the elliptic flow is built up from pressure gradients and 
thus directly proportional to the initial fireball eccentric- 
ity. Unfortunately, there has been no direct experimen- 
tal measurements of this quantity due to the difficulty 
of isolating the initial state contribution from the later 
stages of the fireball evolution. Model estimates of the 
overlap geometry of two nuclei differ up to 20 — 30% in 
eccentricity, which turns out to introduce more than a 
factor of two uncertainty in the extracted values for rj/s 
[r2| . Therefore, the precise determination of rj/s requires 
a more precise knowledge of the initial geometry for the 
produced fireball in the collisions. 

Recently, significant attention has been paid to ini- 
tial geometry fluctuations [2414291 ] which have been used 
to explain the underestimation of elliptic flow calculated 
in various ideal and viscous hydrodynamic simulations 
for the most central collisions. For example, the geome- 
try fluctuations of the positions of nucleons in the Monte 
Carlo Glauber (MCG) model [3(| HH lead to fluctuations 
of the participant plane from one event to another, ren- 
dering larger eccentricities which translates into larger 
elliptic flow for the final state particles. To pursue such 
studies, one needs to run hydrodynamical evolution on 
an event-by-event basis utilizing fluctuating initial con- 
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ditions [32H3J. 

As is known, higher order moments are also present in 
fluctuating initial collision geometry when one performs 
a harmonic/multipole analysis. Triangular geometry and 
flow have recently been proposed to explain features in 
the data such as the ridge and broad away-side corre- 
lations observed in two-particle correlation data, in the 
context of hydrodynamics and transport models I34j36l . 
Higher-order flow coefficients have been measured [371 138[ 
and recent studies show that the initial state density fluc- 
tuations may play an important role in understanding 
the centrality dependence of the ratio v^jv\ [H, To 
achieve a full understanding of the expansion dynamics 
of the produced fireball therefore requires a systematic 
study of initial geometry fluctuations. The main purpose 
of our paper is to investigate how harmonic moments of 
different order propagate through the different stages of 
the fireball history and how they translate themselves 
into the momentum anisotropics of the final produced 
particles. 

In Section II, we construct the full phase space dis- 
tribution of the initial conditions (position and momen- 
tum space) obtained from a Monte Carlo Glauber model 
with the inclusion of the nucleon position fluctuations as 
well as fluctuations from individual nuclcon-nucleon colli- 
sions. The geometry of such initial conditions is analyzed 
in Section III. We study the pre-equilibrium evolution of 
the system and its effect on the spatial geometry in Sec- 
tion IV, where a detailed analysis of the correlations be- 
tween odd and even moments during this period is also 
presented. In Section V, the discretized initial conditions 
is smeared with a Gaussian distribution and, assuming 
sudden thcrmalization, the subsequent evolution of the 
system is modeled utilizing a three-dimensional relativis- 
tic ideal hydrodynamics [3J, [U [42[. Numerical results 
of final state momentum anisotropics after the hydrody- 
namical evolution are presented in Section VI, followed 
by our summary in the last section. 



II. INITIAL CONDITIONS 

Our initial conditions are based on the Monte Carlo 
Glauber model, but differ from other implementations 
of that model as we include the fluctuations of nucleon 
positions as well as the fluctuations originating from in- 
dividual nucleon- nucleon collisions. In addition we ac- 
count for the full phase-space by constructing the par- 
ticle momentum distributions as well. We determine 
the spatial distribution using the well-established two- 
component (binary collision and participant) scaling and 
the momentum distribution is obtained by fitting to data 
on final particle momentum spectra. We also treat the 
early pre-equilibrium expansion of the system using the 
free-streaming approximation prior to the hydrodynami- 
cal evolution. 

We start with the nuclear distribution function inside 



a nucleus taken as the Woods-Saxon form 



Pa{t) 



1 + exp[(r - R)/d] 



(1) 



where the radius R and the diffuse constant d are taken 
as R = 6.38 fm, d = 0.535 fm for a Au nucleus. The 
above distribution is normalized to the atom number 
/ d 3 rp(r) = A with p = 0.163/fm 3 . It is convenient to 
normalize the above distribution function to unity and 
define the single nucleon distribution pa{t) inside a nu- 
cleus, J d 3 rp A (r) = 1. The normalized thickness function 
is defined as 



T A (s) = / dzp A (s,z) 



(2) 



with the normalization J gJ 2 sTa(s) = 1. 

To study the collision between two incoming nuclei at a 
given impact parameter b, one may define the probability 
for a given nucleon i from nucleus A and a given nucleon 
j from nucleus B to collide to be P(si, Sj, b) = cr(sj — Sj — 
b), which is normalized to the nucleon- nucleon inelastic 
cross section <jnn, 



d 2 ser(s) = (7 TV N 



(3) 



where ctnn = 42 mb is taken for nucleon-nucleon col- 
lisions at y/spfN = 200 GeV. From such a probability 
distribution, one may compute the numbers of binary 
nucleon-nucleon collisions and participating nucleons, 



A B 



=J2Y1 [ d 2 s l f A (s l ) f d 2 Sj f B ( Sj )a(s) 
i=i o=i J J 

A 

N palt = £ / ^^Ta^) 

1 - II / dV'flfoOI 1 - ^( s )] j + (A O B) (4) 

where s = Sj — Sj — b. 

To simulate a collision of two nuclei using the Monte 
Carlo approach, one first samples the positions of all nu- 
cleons in the nucleus according to a Woods-Saxon distri- 
bution, and obtains discrete nucleon distributions with 
each single nucleon corresponding to a 6 function. The 
probability function <x(sj — Sj — b) for two nucleons to 
collide is taken to be geometrical in form 



b) = 1 



Si 



b| < y/aNNpK 



a(si - Sj - b) = , |sj - Sj - b| > ^J(t nn /tt (5) 

With this, one returns to the classical picture of col- 
lisions: two nucleons with transverse distance dj_ = 
| Si — Sj — b| < \J<tnn 1^ will collide with each other. 
Note that the assumption of linear trajectories of partic- 
ipant nucleons is still maintained after they collide with 
each other. With the above probability distribution, one 
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reduces the calculation of N co u and N paxi to counting the 
pairs of binary collisions and the number of participating 
nucleons. 

After determining the profiles of two colliding nuclei, 
the produced particle multiplicity in a collisions for a 
given centrality class (or impact parameter b) and ra- 
pidity range A77 can be obtained from the following phe- 
nomenological two-component formula (43[, 



N AA (b,Ari) 



1- 



aJVcou(b) + —-N paxt {b) 



N NN (Ar)) 



(6) 



where Nnn(At]) is the particle multiplicity in a nucleon- 
nucleon collision at the same collision energy. The vari- 
able a controls the balance between two components: 
participant scaling and binary collision scaling. With 
the value of a = 0.13, one may obtain a nice description 
of the centrality dependence of average charged particle 
multiplicity at midrapidity \r)\ < 0.5 in Au+Au collisions 
at ytv]v = 200 GeV (see Fig. Q}. 



• PHOBOS Au+Au 200GeV 
MC-G]aubcr 




FIG. 1: (Color online) Charged particle multiplicity at midra- 
pidity \r]\ < 0.5 as a function of centrality in Au+Au collisions 
at y^W = 200 GeV. 

As is well known, the particle multiplicity in high en- 
ergy collisions is fluctuating from one event to another. 
The distribution of particle multiplicities N(Arj) for a 
given rapidity range A77 may be well described by a neg- 
ative binomial (NB) distribution, (HrH?]], 



• UA5 p+p 

STAR p+p 

♦ NB distribution 



FIG. 2: (Color online) Charged particle multiplicity distri- 
bution at midrapidity \r]\ < 0.5 in p+p collisions and p+p 
collisions at ^/snn = 200 GeV. 



For high energy nucleus-nucleus collisions, we include 
particle multiplicity fluctuations by evaluating Eq.Q 
on an cvent-by-event basis, with the distribution of 
Nnn(At]) given by Eq.((7]). The balance factor a of bi- 
nary collision and participant scaling in Eq.(j6|) is imple- 
mented by randomly keeping only a fraction a of particles 
from a binary collision, and a fraction (1 — a)/2 of parti- 
cles originating from a participating nucleon. In general 
the positions of particles might be sampled according to a 
smeared distribution around the positions of binary colli- 
sions or participant nucleons. Here we take the positions 
of particles as the same positions as binary collisions or 
participant nucleons. The Lorentz contraction in the lon- 
gitudinal direction is taken into account by contracting 
the longitudinal position z of each particle by a factor of 
100 for Au+Au collisions at ^/snn = 200 GeV. 

With the number of produced particles and their posi- 
tions fixed, we also assign momenta to each particle. In 
this work, particle transverse momenta px are sampled 
according to the following power law distribution, 



dN 
dpi, 



(8) 



P(N,fi,k) = 



T{N + k) {n/k) N 
T(N + l)T(k) (n/k + l) N + k 



(7) 



where /1 is the mean of the distribution, and k is related 
to the shape of the distribution. The variance is given 
by a 2 = n(fJ,/k + 1) and the scaled invariance is defined 
as uj = a 2 /fi = fi/k + 1. With the values of /1 = 2.35 and 
k = 1.9, one obtains a good description of the charge 
particle multiplicity measurements for both p+p colli- 
sions from UA5 [45[, and p+p collisions from STAR [46| 
at midrapidity \r/\ < 0.5 (see Fig. [5]). 



where a is the normalization constant, and b and c are 
taken as b = 0.88 and c = 4. The azimuthal angle of the 
transverse momentum is uniformly distributed. Particle 
rapidities are taken to be uniformly distributed around 
mid-rapidity \y\ < 1 (When changing the rapidity range, 
we change both the mean /1 and the parameter k and keep 
the scaled variance w of NB distribution fixed). Particles 
with large rapidities will be absent from the central ra- 
pidity region at later times and we neglect them in this 
work. With the above setup, we obtain the full phase 
distribution of the system at initial production time. 
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III. INITIAL GEOMETRY 

Before moving to the evolution of the system, we first 
investigate its geometrical properties at the production 
time. In a nucleus-nucleus collision, the reaction plane is 
defined by the beam direction (z) and the impact param- 
eter direction (x). The impact parameter direction and 
the third orthogonal direction (y) define the transverse 
plane (one typical collision event is shown Fig. [3]). We 
call the plane defined by z direction and y direction the 
vertical plane. The geometry of the transverse plane is 
particularly interesting due to the fact that the elliptic 
flow V2 is found in ideal hydrodynamics to be propor- 
tional to the initial eccentricity ti of the overlap region 
of the colliding nuclei, the determination of which plays 
an important role in the extraction of the transport co- 
efficients of the produced fireball. 

For averaged initial conditions, the geometry of the 
system can be studied directly in the above framework 
due to the coincidence of the vertical plane and the spa- 
tial event plane for e2 (a rotation by 7r/2 of the partici- 
pant plane if the participating nuclcons are considered for 
the spatial distribution). With fluctuating initial condi- 
tions, the spatial event plane is tilted with respect to the 
reaction plane from one event to another. We call the 
angle between the spatial event plane and the reaction 
plane the spatial event plane angle <f>2- Note that this 
choice of the spatial event plane is convenient when gen- 
eralizing to higher moments since the event plane angle 
distribution always has a maximum in the y direction for 
all even moments whereas not always one of the minima 
is in the x direction. 

The final elliptic flow vi is defined with respect to a 
third plane, the momentum event plane that is recon- 
structed in experiments from the measured momentum 
distribution of the produced particles. Again in ideal 
hydrodynamics with smooth initial conditions this event 
plane coincides with the reaction plane, i.e., is rotated 
with respect to the spatial event plane by 7r/2. This ro- 
tation ensures that the final vi has the same sign as the 
initial 62. If an event-by-event analysis with fluctuating 
initial conditions is applied, a strong correlation of the 
final momentum event plane to the initial spatial event 
plane still remains, but fluctuates around n/2 as has been 



shown in [33 1 



One may generalize the above concept for every har- 
monic moment, and define the spatial anisotropy pa- 
rameters e„ as follows. The first moment t\ can al- 
ways be made to vanish by shifting the coordinates to 
the center of mass (CM) frame of the system such that 
(x) = (y) = 0. Note (...) throughout this paper represent 
averages over the phase space profile for a given event, 
except in the Appendix. Once the system is shifted to 
its CM frame, t\ = 0, all higher harmonic moments can 
be defined as: 




FIG. 3: (Color online) The transverse plane for one typi- 
cal collision event, where the cycles represent nucleons from 
two nuclei, with shaded ones for participating nucleons. Also 
shown are the locations of different planes: the reaction plane 
(RP), the spatial event plane (SEP) and the momentum event 
plane (MEP) for n = 2. 



where r± = y 1 x 1 + y 2 , and tfi = arctan(y/a;) are polar 
coordinates in the transverse plane. The spatial event 
plane angle with respect to the reaction plane can be 
found through the following formula, 



4\ 



I (r'f am(nd>)) 

- arctan ±+ )— 

n (r" cos(n<p)) 



(10) 



Note in our definition $„ fluctuates from one event to an- 
other in the range of (— ir/n, 7r/n), but it is equivalent to 
rotate such angle by 2n/n. Once the event plane angle is 
found, the definition of the spatial anisotropy parameters 
may be reduced to e„ = (r™ eos[n.(<^> — )])/(?"")• 

The flow coefficients v n are defined as the n-th Fourier 
moment of the particle momentum distribution with re- 
spect to each momentum event plane, 



v n = (cos[n(V> - *„)]) 



(11) 



yVl cos(n0)) 2 + (rj sin(n</>)) 2 / (r\) 



(9) 



where ip = tan -1 (p y /p x ) is the azimuthal angle of par- 
ticle momentum p in the CM frame. Here we define the 
momentum event plane by a rotation angle n/n with re- 
spect to the initial spatial event plane, ^ n = + n/n. 
This rotation is just a convention generalized from the re- 
quirement that a positive initial eccentricity generates a 
positive value of final elliptic flow in ideal hydrodynamics 
with average initial conditions. Note our definition of the 
momentum event plane does not necessarily corresponds 
to the event plane that is reconstructed in experiments as 
just mentioned, but for our systematic study it provides 
an unambiguous basis to quantify the final state response 
to the initial state anisotropies. 

Since the event plane defined by Eq. (fl0|) fluctuates 
around the vertical plane (y direction) for even moments, 
we may perform a transformation <f>„ — > <$>' n by evaluating 
Eq. (flUl) with — >• (ft' = (p + 7T /2. This corresponds to a 
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*n" 



FIG. 4: (Color online) The probability distribution of $^ at 
production time with b = 8 fm. 



rotation of the coordinate system (x, y) by 7r/2 which en- 
sures the distributions of all even $^ peak at as shown 
in Fig. 2J In this plot, the impact parameter b is taken 
to be 8 fm for all events. While all even moments are 
strongly correlated with the reaction plane with one max- 
imum along y direction, all odd moments arc uniformly 
distributed. This may be understood since the odd mo- 
ments of spatial anisotropy purely originate from fluctu- 
ations while the even ones are combined effects of fluc- 
tuations and geometry. As a consequence, if one defines 
the spatial anisotropy parameters e„ with respect to the 
pre-determined the reaction plane, the event-averaged e„ 
vanishes for all odd moments, but not for even ones. We 
also observe that the distributions of even moments is 
wider for higher values of n due to weaker correlations 
with respect to the reaction plane [in fact, what mat- 
ters is the distribution of n$^, as $^ fluctuates within 
(— ir/n, tt/ti)]. 



tributions are plotted as a function of impact parameter 
b. One can see that the widths all odd values of n align 
with each other at as expected [for a uniform distri- 

bution from —ir/n to tv/ti, the variance is a 2 = ir 2 /(in 2 ) 
and we are plotting na(Q' n ) = a(n$>' n )]. Also due to 
symmetry in central collisions there is no correlation be- 
tween the angle <f>^ and y direction for all values of n. The 
anisotropy is purely from fluctuations, rendering uniform 
distributions also for even values of n. As one moves to 
non-central collisions, geometry comes into play and may 
dominate over pure fluctuations, hence the widths of even 
n distributions become smaller. For very peripheral col- 
lisions, the importance of the geometry diminishes due to 
the small size of the system, and even n distributions be- 
come broader again. We also observe that even harmonic 
moments with higher values of n have weaker dependence 
on centrality. 




b (fm) 



FIG. 6: (Color online) The spatial anisotropy parameters e n 
at production time as a function of centrality from Monte 
Carlo Glauber modelling of Au+Au collisions at ^/snn = 
200 GeV (see text for details). 
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FIG. 5: (Color online) The widths of <&' n distribution (times 
n) at production time as a function of impact parameter b. 

We also investigate the centrality dependence of the 
above correlations in Fig. [5J where the widths of the dis- 



In Fig. [51 the first few spatial anisotropy parameters e n 
are plotted as a function of impact parameter b (ei = 
is not shown). One may observe that all moments are 
of the same magnitude for typical non-central collisions. 
In central collisions, as pure fluctuations instead of ge- 
ometry generate the anisotropy, higher moments acquire 
larger values due to larger fluctuations brought by the 
power n in the definition of e„. Note that if the same 



weight, i.e. 



is taken for every e„ as in Ref. [3£ 



all moments are the same in central collisions and £2 is 
larger than all higher moments in non-central collisions 
(also note e\ is non-zero if r 2 ^ is used). 

In our initial conditions, we may also calculate the mo- 
mentum anisotropics as we have generated the full phase 
space distribution for the produced system. Since the 
initial particles are sampled with a symmetric azimuthal 
distribution, one obtains zero v n when averaging over 
events. In Fig. [7J the width of the initial v n distribution 
is plotted as a function of impact parameter b. We find 
that the width increases as one moves from central colli- 
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FIG. 7: (Color online) The widths of the initial v„ distribution 
at the production time as a function of impact parameter b. 



sions to non-central collisions due to the decrease in the 
number of particles in the produced system. In fact, the 
v n distribution is a Gaussian as a result of central limit 
theorem and its width is found to be 1/ %/2A for all val- 
ues of n, except for very peripheral collisions where the 
particle number N is too small. The figure contains two 
sets of curves: the upper one is for all particles within a 
rapidity bin of \y\ < 0.5 and the lower one for \y\ < 1. 
The curves are related by a factor of \/2 since the num- 
ber of particles in the system is doubled via the doubling 
of the rapidity bin. The non-zero width of initial v n dis- 
tribution helps to explain the wide distribution of the 
transformation matrix elements between initial e„ and 
final v n as shown in Fig. [TSJ It serves as another source 
that contributes to final flow fluctuations in addition to 
initial state geometry fluctuations. 



IV. P RE-EQUILIBRIUM PHASE 

As to date, hydrodynamical simulations mostly use ini- 
tial conditions calculated at the initial production time 
of the medium, and have neglected the influence of the 
pre-equlibrium time evolution of the colliding matter. In 
this sense, the spatial information inferred from com- 
paring experimental measurements with hydrodynamical 
simulations is for the system at the starting time of the 
hydrodynamic evolution t = to, not at the initial produc- 
tion time t = 0. However, the pre-equilibrium evolution 
may be important to include when considering the geom- 
etry fluctuations of the produced matter. Unlike the hy- 
drodynamical evolution which directly translates the ini- 
tial geometric anisotropics into the observed momentum 
anisotropics, the early pre-equilibrium expansion of the 
system will not only smear out the spatial fluctuations 
and change the local momentum distribution, but may 
also lead to correlations between odd and even moments. 
The inclusion of the pre-equilibrium evolution could be 
also important for studying the Hanbury-Brown-Twiss 



intcrfcrometri radii as it will generate some amount of 
early flow gHHl- 

To simulate the pre-equilibrium evolution, we solve 
the Boltzmann equation for the phase space distribution 
/ (x, p, t) = dN/d 3 xd 3 p of the system, 



(d t + vV x )/(x,p,t) = C[/] 



(12) 



Here, for simplicity massless particles are considered, 
|v| = 1. The free-streaming term will smear out the 
spatial fluctuations and change the local momentum dis- 
tribution due to pure fluctuations. The collision term is 
important for studying the details of the thermalization 
of the system, a complex issue which has not been fully 
understood yet. In this work, we focus on the effect of the 
pre-equilibrium expansion on the system geometry and 
only include the free-streaming term by setting C[f] = 0. 
Such a treatment is important for our study as the flow 
built up by hydrodynamical evolution is mostly driven 
by spatial anisotropy. The consideration of the collision 
term will remain for a future project. The Boltzmann 
equation containing only the free-stcaming term can be 
solved analytically, with the streaming solution given by 
/(x,p,i) = /(x — v(t — h),p, U), where t{ is the initial 
starting time of the evolution. 



0.95i - ■ 

-i ■ 





— # n 


= 2 




-■ n 


= 3 


♦ - 


■♦n 


= 4 


tr- 


■in 


= 5 


4 


M n 


= 6 




6 

b (fm) 



FIG. 8: (Color online) The ratio of e n at to = 0.6 fm/c to those 
at initial production time as a function of impact parameter 
b. 

The effect of free-streaming on the spatial anisotropics 
during the early expansion is shown in Fig. [8j where 
the ratios of anisotropy parameters e„ evaluated at to = 
0.6 fm/c to those at the production time are shown as a 
function of centrality. As expected, the expansion of the 
matter due to free-streaming smears out the spatial fluc- 
tuations: all the spatial anisotropy parameters e„ become 
smaller. This diminishing effect is more pronounced in 
non-central collisions due to the smaller size of the sys- 
tem. We also observe that higher moments get more 
diminished than lower moments. 

We further explore the time evolution of the above 
smearing effect due to free-streaming as shown in Fig. 
[21 where e n at an impact parameter of b = 8 fm is 
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FIG. 9: (Color online) Time evolution of e n when particles are 
just freely streaming for an impact parameter of b = 8 fm. 



plotted function of time. We observe that the 

anisotropy parameters e ra decrease rather fast for the first 
2 — 3 fm/c and then slowly saturate. The observed re- 
duction of the e„ hints at the importance of including 
the pre-equilibrium expansion when studying the initial 
state geometry fluctuations for hydrodynamical simula- 
tions. The relative size of the different coefficients even 
depends on the duration of the pre-equilibrium expan- 
sion and the relative size of the resulting flow coefficients 
might be used to constrain this initial evolution. 

As we just mentioned, there are two separate effects 
due to the pre-equilibrium evolution: the pure drift effect 
due to the system expansion which tends to diminish all 
moments, and the correlations between odd and event 
moments which is absent in the later hydrodynamical 
evolution. The mixing effect between odd and event mo- 
ments can be clearly seen when one explicitly performs 
a multipole-expansion analysis for the Boltzmann equa- 
tion. To start, we choose spherical polar coordinates for 
both coordinate space x = (r, 9, <j)) and momentum space 
p = (p,O p ,(j>p). To make the analysis dimcnsionless, we 
define f = r/r max and p — p/T p , where r max and T p are 
two dimensional quantities (being constants or varying 
with time). Then we expand the phase space distribu- 
tion as 



/(x,p,t)= ]T a^J 1 {t)R nl {a nh f)Y lmK 

nlmNLM 

eM-p)PN(p)Y L u{0 P Ap) (13) 
In the above expression, 



ji+i(a n i) 



Pn(P) 



(iV + /x)! 



(14) 



where ji is the spherical Bcssel function with a n i the nth 
root of function ji , the /i-th order Laguerre function 



(Here we choose /j = 2) and Y" l (8, <fi) the spherical har- 
monics. The expansion coefficients a^h^ 1 are determined 
from the phase space distribution by 



a 



NLM, 



pi p pOC p 

i (*) = / r 2 dr I dVL I p 2 dp / dfl p R n i(a n i,f) 

r^(fi)p JV (p)y J j M (n l ,)/(x,p,t)(i5) 



The spatial anisotropy parameters e m are related to the 

NLM 
nlm 

C\m] 



expansion coefficients a^^ M , by 



((sin0) m sin(m<7i)) 

((sin(9) m cos(m0)) 
where 

cm = (-i) 

J r [n, m 



- L J nmmJ Jr[n, TTl] 



C[m 



N 



E Im K 



!X>K™J MnM (16) 



47r(2m)! V8tt 
2m + 1 (2m- 1)!! 

V2 



f drj m (a nm f) (17) 



The normalization factor N represents the total number 
of particles in the system, which is given by 



N 



(18) 



Note for ((r sin^)" 1 cos(to0)) and ((r sin6>) m sm(m(f>)) 
which appear in the definition of e m , there will be an 
extra factor r m in the integral J r . 

Within the above multipole expansion analysis, the 
Boltzmann equation becomes, 



in 



E a^ m M (t)^Ll r [n',l',n,l] 

nlmNLM max 

Snn' (8i,i>-i + Si i'+i) (5l,l>-i + 5 L ,L'+\) 



l + l' + l 2L + 1 



-(L,0;l,0|i',0) 



y 2(2/ + 1) V 2L' + 1 
y~*d m ,m'+iSM,M'-i(l', m'; 1, i\l, m)(L, M; 1, i\L', M') 

i 

= C[a»;^'\ (19) 

where 

V2 V2 



7 r [n', n, 1} = 



jl'+i(a n 'i') ji+i(a n i) 



f 2 dfj v (a n >i>r)jv {a„if) 



(20) 



Note the indices in ji>(a n ir) in ijP^ [n',l',n, I], which do 
not allow us to perform the integral using orthogonal re- 
lations. (Ji, rm; 1-2, m 2 |j, m) are the Clebsch-Gordan coef- 
ficients for adding two angular momenta j = lx + 1 2 - More 
details of the derivation are presented in the Appendix. 
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FIG. 10: (Color online) Time evolution of ((sin 6) n cos(n0)} 
from multipole expansion analysis (lines) and from directly 
solving the free-streaming term of Boltzmann equation (sym- 
bols) for an event with 6 = 8 fm. 




FIG. 11: (Color online) The probability distributions of the 
transformation matrix between £2 and £3 with the early ex- 
pansion time taken as to = 1.2 fm. The numbers represent 
the mean of each distribution. 



From the above equations, one immediately sees the 
mixing between odd and even moments both for spatial 
part (I — > I ± 1) and momentum part [L — > L ± 1) due to 
the free-streaming of particles. We also check the above 
expression by comparing numerically with the result from 
directly solving free-streaming part of Boltzmann equa- 
tion. This is shown in Fig. [TU1 where we plot the time 
evolution of ((sin#) n cos(n0)) for one typical event with 
impact parameter b = 8 fm and see that the two results 
nicely agree with each other. 

In order to further separate the correlation effect from 
the pure drifting effect during the pre-equilibrium expan- 
sion, we perform the following analysis. We relate the 
spatial anisotropics at two different times by a transfor- 
mation matrix, 



C2(to) 

e 3 (*o) 



#22(^0) 

D 32 (t ) 



D 2S (t ) 
D 33 (t ) 



e 2 (0) 
ca(0) 



(21) 



Here we only consider the second and third moments - 
the inclusion of higher order moments is straightforward 
and is expected to give only small contributions which we 
neglected in the current analysis. The diagonal elements 
of the transformation matrix quantify the pure drifting 
effect and the off-diagonal elements represent the effect 
of the mixing between the second and third moments. 
We obtain the distribution of the transformation matrix 
elements by pairing two linear independent events from 
a large set of events. 

In Fig. II H we show the probability distribution of four 
elements of the above transformation matrix, with the 
pre-equilibrium expansion time taken to be to = 1.2 fm/c. 
The numbers in the figure represent the mean of each dis- 
tribution (to which the arrows point). The impact pa- 
rameter is taken to be 8 fm for all events in this plot. One 
clearly observes the smearing effect from the free stream- 
ing when one looks at the distributions of the two diago- 
nal elements. The pure drifting effect is more pronounced 



for the third anisotropy parameter 63 (19%) than for the 
second one €2 (13%), consistent with the above results. 
The two off-diagonal elements are close to zero imply- 
ing weak correlations between e 2 and £3 originating from 
the free-streaming of the system. We further investigate 
the time evolution of these matrix elements up to 2 fm/c 
in Fig. [T2] Both the drifting effect and the mixing of 
even and odd moments tend to increase with time as the 
system expands. 
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FIG. 12: (Color online) Time evolution of the transformation 
matrix between £2 and £3. 



V. HYDRO DYNAMICAL EVOLUTION 

In the previous sections, we have presented the initial 
conditions of the system at production time and simu- 
lated the pre-equilibrium evolution by utilizing the free- 
streaming approximation. As we have not included in- 
teraction among the produced particles, the system is 
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still highly non-thermal. Up to know, little knowledge 
has been attained about the details of the thermaliza- 
tion mechanisms in relativistic heavy-ion collisions. In 
this work, we follow the common practice to assume a 
sudden thermalization of the system at t = to and start 
the hydrodynamical evolution with the initial conditions 
obtained above (including the free streaming evolution 
from t = to t — to). We first calculate the energy- 
momentum tensor from the full phase space distribution 
/(x,p,i), 



T" v {x) 



J ^y7(x,p, 



t) 



(22) 



For our discretized phase space distribution /(x, p,t) = 
J2 i 6(x — Xi)J(p — Pi), the momentum integration J d 3 p 
turns into sums over all particles. The discretized spatial 
part is smeared with a Gaussian function in order to en- 
sure a sufficiently continuous distribution necessary for 
the hydrodynamic simulation, 



<5(x - Xi) 





(x-Xi) 2 +(j/-3/i) 2 " 




[ (*-*«)' 1 


exp 


exp 


3*3, 


2ai \ 



(23) 



where the widths a xy and a z characterize the granularity 
of the system in the transverse and longitudinal direc- 
tions. Physically, this procedure can be interpreted as 
thermal smearing of the system which should have oc- 
curred prior to thermalization at t = to- In general the 
choice of these width parameters depends on the dura- 
tion of the pre-equilibrium phase and the thermalization 
time at which one starts the hydrodynamic evolution. 
Different choices of the smearing width will affect the lo- 
cal density of the system, and thus influence the spatial 
anisotropy parameters. 




0.4 0.6 
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FIG. 13: (Color online) The ratio of e n with smearing to those 
without smearing as a function of the transverse Gaussian 
width. 

The effect of the Gaussian smearing on the spa- 
tial anisotropy is shown in Fig Q21 where the ratio of 
anisotropy parameters e„ with smearing to these without 
smearing is shown as a function of transverse smearing 



width cr xy . As we are studying the spatial anisotropy 
in the transverse plane, the smearing of the longitudi- 
nal direction should be irrelevant and we fix it to be 
o z = 0.5 fm for our study. The impact parameter is 
taken to be 8 fm for all calculations shown in this fig- 
ure. As expected, the spatial anisotropics are reduced 
as one increases the width of the transverse Gaussian 
function. Similar to pre-equilibrium evolution shown be- 
fore, such smearing effect is more prominent for higher 
moments than for lower moments. Combining both ef- 
fects (pre-equilibrium evolution and Gaussian smearing), 
for typical non-central collisions t2 may be reduced by 
about 10% for a Gaussian width of a xy = 0.5 fm and a 
typical pre-equilibrium evolution time of to = 0.6 fm; a 
factor of 2 larger effect is observed for for £4. 

In the above construction of the energy-momentum 
tensor, the initial conditions at production time are fit- 
ted to the final state particle multiplicity distribution at 
midrapidity (see Fig. [T]and[2]). Therefore, the energy 
density of the system at the thermalization time to is 
underestimated, due to the longitudinal (and transverse) 
expansion during the hydrodynamical evolution. This 
effect can be estimated to be about a factor of 2.2 by 
directly comparing our calculation to the final average 
charged particle multiplicity dN^/dr) w 700 in central 
Au+Au collisions at ^/snn = 200 GeV. We have not 
tuned our parameters to match the final state particle 
spectra as we are here not aiming at providing a com- 
prehensive quantitative description of the time-evolution 
of a heavy-ion collision, but rather at a targeted study 
of initial state fluctuations and how these initial spatial 
anisotropics propagate through the fireball history and 
translate themselves into collective flow in the final state. 

In Fig. Q3J we show the a few snapshots of the energy 
momentum tensor components in the transverse plane 
(the horizontal and vertical axes are x and y axes in unit 
of fm) for one typical event with an impact parameter 
b = 8 fm. To illustrate the effects of the pre-equilibrium 
evolution and Gaussian smearing of discretized space dis- 
tribution, we plot three different sets of pre-equilibrium 
time and Gaussian smearing width: the upper for to = 
and a xy = 0.5 fm, the middle for to = 0.6 fm/c and 
o'xy = 0.5 fm, and the lower for to = 0.6 fm/c and 
a xy = 1 fm. On the left we show the distribution of T 00 
and on the right the flow vector (T 0x /T 00 , T°« /T 00 ), with 
the arrows representing the directions and the lengths of 
arrows for the relative magnitudes of the vectors (within 
each plot). Comparing the left and middle panels, one 
can clearly see that pre-equilibrium evolution makes the 
system larger (thus the energy density becomes smaller) 
and generates some amount of radial flow. The effect of 
Gaussian smearing can be seen by comparing the middle 
and right panels: both the energy density and the flow 
velocity smoothen out significantly when one increases 
the Gaussian width. 

After obtaining the energy momentum tensor as de- 
scribed above, we directly start the hydrodynamical evo- 
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FIG. 14: (Color online) The distributions of the energy mo- 
mentum tensor components in the transverse plane (x, y) for 
one typical event with impact parameter b = 8 fm: the left 
for T 00 and the right for the flow vector (T 0x /T 00 , T 0y /T 00 ). 
Three different sets of parameters are used: to = 0, a xy = 
0.5 fm (upper), to = 0.6 fm/c, a xy = 0.5 fm (middle), and 
to — 0.6 fm/c, G xy = 1 fm (lower). 



lution with the assumption of a sudden thcrmalization, 



d^ v {x) = o 



(24) 



Here an ideal hydrodynamical evolution code [4l|, |42( is 
utilized for our study with a lattice equation of state 
[52l [53| for the hot and dense matter created in Au+Au 
collision at tJsnn = 200 GeV. Particle production at the 
end of the hydrodynamic evolution when the matter is di- 
luted in the late stage is treated as a gradual freeze-out on 
an approximated iso-eigentime hyper-surface according 
to the Cooper- Frye prescription y52|, [f34| . For simplicity, 
we have not taken into account the hadronic rescattering 
in the dilute hadron gas and the resonance decays, since 
they should not have much influence on the results for 
the charged particle flow coefficients as has been shown 
in 13211 . 



VI. FROM INITIAL GEOMETRY 
FLUCTUATIONS TO FINAL FLOW 

The above event- by- event setup of the system evolu- 
tion from initial production time to the freeze-out of the 
final state should include all ingredients that are neces- 
sary for the study of the build-up of collective flow during 
the hydrodynamical evolution. For the following results, 
we use the produced charged particles with transverse 
momenta pt < 2 GeV/c and peudorapidity |?y| < 1. 
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FIG. 15: (Color online) The first few spatial anisotropiy pa- 
rameters e n and flow coefficients v n as a function of n for 
b = 5 - 10 fm. 

In Fig. [TSl we show the first few flow coefficients v n for 
final state particles together with three different spatial 
anisotropics evaluated at the production time, at t = to 
before and after Gaussian smearing. In this figure, the 
impact parameter is randomly sampled in the 5 — 10 fm 
bin according to the probability distribution P(b) oc 6, 
the pre-equilibrium evolution time is set as to = 0.6 fm/c 
prior to the hydrodynamical evolution, and the Gaussian 
widths for smearing the discretized initial conditions are 
taken to be a xy = a z = 0.5 fm. One can see that all 
flow coefficients v n with n greater than 5 are negligible 
and only the first few v n (n = 2, 3, 4) survive after the 
hydrodynamical evolution. We may conclude that the 
analysis of flow coefficients v n allows only for the extrac- 
tion of the first few spatial anisotropy parameters e n , 
but may not provide sufficient information to recover the 
full initial geometry in terms of all of its higher order 
harmonics. In order to achieve this, one needs additional 
observables with an increased sensitivity to the higher or- 
der spatial anisotropy parameters. We also note that the 
spatial asymmetry obtained from comparing flow mea- 
surements to a hydrodynamical simulation only applies 
to the spatial characteristics of matter at the starting 
time of hydrodynamical evolution ty. In order to obtain 
the geometry at the initial production time t = 0, one 
needs to account for the dynamics of the pre-equilibrium 
phase - in our analysis this would be the smearing effects 
during the initial free streaming of the particles and the 
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smoothing of the discrctizcd initial conditions as shown 
in the figure. 




0.1 



FIG. 16: (Color online) The first few flow coefficients v n as 
a function of the same order spatial anisotropy parameter e„ 
for b = 5 - 10 fm. 



In order to study the response of flow build-up to 
the initial geometry, Fig. [16] shows the first few flow 
coefficients function of the corresponding spa- 

tial anisotropy parameter e„ evaluated at the production 
time. As expected, the hydrodynamical evolution trans- 
lates spatial anisotropics into momentum anisotropics, 
resulting in an essentially linear relation between v n and 
e n . We also observe that the curves have smaller slopes 
for higher moments and become more or less fiat when n 
is greater than 4 — 5. This implies that higher flow coef- 
ficients show weaker response to the corresponding geo- 
metrical harmonic moments due to larger diminishing ef- 
fect originating from the combination of pre-equilibrium 
evolution, Gaussian smearing of the discretized spatial 
distribution, and the hydrodynamical evolution, all of 
which tend to lead to a larger suppression for the higher 
moments. 
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FIG. 17: (Color online) The first few flow coefficients v n as a 
function of the second order spatial anisotropy parameter t2 
for 6 = 5 - 10 fm. 



We also explore the correlations between different har- 
monics moments. As an illustration, we plot the first few 
flow coefficients function of the second spatial 

anisotropy parameter ti. On can see that except for i>2 
versus £2, all other curves are essentially flat due to the 
combined effect of the small correlations between odd and 
even moments and the reduced effect on higher moments 
from the pre-equilibrium and hydrodynamical evolution. 





-• P(M 22 ) 


■- 


-■ P(M,J 


♦ - 


■♦ P(M 25 ) 




■A P(M„) 




FIG. 18: (Color online) The probablity distributions of the 
transformation matrix between the initial spatial anisotropy 
parameter £2,3 and the final flow coefficients 112,3 

In order to separate the pure smearing effect from the 
mixing effect between different moments, we perform an 
analysis similar to that for the pre-equilibrium evolution. 
On may define the transformation matrix M nm between 
the initial spatial anisotropy and the final flow coefficients 
as 



M 



nm u m 



(25) 



where M nm characterizes the strength of the coupling 
between initial e m and final v n . The diagonal elements 
of the transformation matrix quantify the response of v n 
to e„ and the off-diagonal elements represent the effect of 
the mixing response between different moments v n and 
e rn . Here again we only consider the second and third 
moment, 



M22 M 23 \f e 2 
M 32 M 33 M e 3 



(26) 



The extension of this ansatz to include higher order mo- 
ments is straightforward and expected to only give small 
contributions to the dominant moments n = 2,3. The 
distribution of the transformation matrix elements is ob- 
tained by solving two linear independent equations which 
correspond to a pair of linear independent events chosen 
from a large set of events. 

In Fig. [TH we show the probability distribution of the 
four elements of the transformation matrix between final 
V2, V3 and initial £2, £3- We find that for the two diagonal 
elements (A^22)cvt = 0.21 is larger than (-M33) ev t = 0.13, 
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implying stronger response of vi to £2 than ^3 to £3 as 
expected from Fig. [T5] The two off-diagonal elements 
(-A^32)cvt and (M23) cv t again are very small, implying a 
weak response of ^3(^2) to £2(63) during the hydrodynam- 
ical evolution. Another interesting feature is the wide 
distribution of the transformation matrix which encodes 
the fluctuations of final flow coefficients. We note two 
initial state effects that contribute to such wide distri- 
bution: initial geometry fluctuations and initial v n fluc- 
tuations (see Fig. [7]). If one wants to extract informa- 
tion about the initial collision geometry from the mea- 
sured flow anisotropics, it is important to separate the 
two sources of fluctuations in the transformation matrix. 
Experimentally, this could be achieved by measuring the 
rapidity correlations of the final flow coefficients since 
the initial state geometry fluctuations are expected to 
be long-range in rapidity, while initial state flow fluc- 
tuations should decrease when one makes the rapidity 
window wider. 
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FIG. 19: (Color online) The final flow coefficients 
function of impact parameter b. 



Finally, we explore the centrality dependence of the fi- 
nal state flow coefficients as shown in Fig. [19j The split- 
ting of v n is clearly seen for all centralities: the lower 
v n are larger than higher v n , at least for the first few v n 
(when n > 5, v n are so small that it is difficult to re- 
solve their splitting). For the most central collisions, our 
statistics is not sufficient to distinguish between different 
curves, but one would expect the same ordering, even 
though the splitting might be smaller. This is due to pure 
fluctuations being the only source of spatial anisotropics 
for both odd and event moments in central collisions, 
and the smearing effect in the subsequent pre-equilibrium 
evolution is more prominent for higher moments, thus 
leading to less flow builtup for higher order v n during 
the hydrodynamical evolution. One observes different 
centrality dependencies for different v n : the lowest v n 
have the strongest centrality dependence. 



VII. SUMMARY 

We have presented a systematic study of initial col- 
lision geometry fluctuations and have investigated how 
they evolve throughout the whole history of the colli- 
sion and finally translate into measurable momentum 
anisotropics of the produced hadrons. Our initial condi- 
tions at production time t = are obtained via a Monte 
Carlo Glauber model with the inclusion of nucleon posi- 
tion fluctuations, plus additional fluctuations stemming 
from individual nucleon-nucleon collisions. In addition 
we make an ansatz for the initial transverse momentum 
distribution of the produced particles which is impor- 
tant for the treatment of the pre-equilibrium phase of 
the collision. We evolve our full phase space distribu- 
tion using a Boltzmann equation and approximate the 
pre-equilibrium evolution by treating all particles as free 
streaming. A sudden thermalization of our initial condi- 
tions is enforced for the subsequent hydrodynamic evolu- 
tion of the thermalized system, which is performed using 
three-dimensional relativistic ideal hydrodynamics. 

Our analysis shows that though all initial spatial 
anisotropy parameters are of the same magnitude, 
only the first few flow coefficients for the momentum 
anisotropy of final state hadrons actually survive after 
hydrodynamical evolution. We also quantitatively inves- 
tigate the mixing between odd and even harmonic mo- 
ments during the pre-equilibrium evolution and its ef- 
fect on the evolution of the system asymmetry is found 
to be small. The anisotropy of the matter is found to 
be affected by the pre-equilibrium evolution and by the 
smoothing of the discretized initial conditions necessary 
for the hydrodynamical evolution (which can be seen as 
equivalent to thermal smearing expected to occur during 
thermalization) , both of which tend to smear out the spa- 
tial anisotropics. The hydrodynamical evolution leads to 
an additional dampening of the flow response to the ini- 
tial spatial anisotropics, particularly for the higher order 
moments. This makes it difficult to recover the full ge- 
ometry of the initial state from measuring high order flow 
coefficients which have the ability to provide additional 
information on the transport properties of the produced 
matter. We also observe the contribution of initial state 
flow fluctuations to final flow fluctuations, which could 
be separated by rapidity correlation measurements for 
a better understanding of initial state geometry fluctua- 
tions. 

In summary, we have conducted an event-by-cvent 
study of the time evolution of the multipolc moments 
of the collision geometry in a relativistic heavy-ion col- 
lision. Our study sheds light on how these multipolc 
moments relate to measurable collective flow coefficients 
of the hadronic final state and how the collision dynam- 
ics, both in the pre-equilibrium and in the hydrodynamic 
evolution phase, affect the correlation between the ini- 
tial spatial anisotropics and the final momentum space 
anisotropics. It allows for improved constraints on the 
determination of various transport properties of the QCD 
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medium, which commonly are extracted by analyzing the 
observed momentum anisotropy of the final particles and 
are very sensitive to the proper description of initial spa- 
tial anisotropics. Our current work can be improved in 
many directions. Here, we have only focused on the qual- 
itative study of the propagation of the initial state geom- 
etry fluctuations; a more quantitative study including a 
comparison with experimental measurements would be 
desirable. The pre-equilibrium phase has been approxi- 
mated by free-streaming of particles; we fully expect that 
the inclusion of a realistic collision term in the Boltz- 
mann equation would provide more sophisticated initial 
conditions for the hydrodynamic evolution. We have 
performed our calculations using ideal hydrodynamics; 
improving these with the use of viscous hydrodynamics 
should help to separate viscosity dominated effects from 
non- viscous effects on the evolution of the geometry fluc- 
tuations. All these tasks we leave to future investigations. 



The second term involves the gradient of a function of 
r times a spherical harmonics. We note the following 
gradient formula, 



VF(r)Y lm (f) = J2ii 



I f dF l + l 



21 + 1 V dr 



-F 



Yi- ltm -i(f)(l -l,m-i, l,i\l,m) 
dF I 



l + l 
21 + 1 



F Yi + i tTn -i(r)(l + 1,to - i, l,i\l,m) 



(A4) 



Here (Zi, mi, I2, m.2 1 J, m) represent Clebsch-Gordan coef- 
ficients for adding two angular momenta j = li + I2. The 
spherical basis vectors are defined as 



£±1 = +^(ex±ie y ), io = e z 



(A5) 
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Appendix A: Multipole Expansion 

In this appendix, we present more details of the mul- 
tipole expansion of the left hand side of the Boltzmann 
equation. The first term is straightforward, 



at 



= E 



nlmNLM 



Of 



Rnl(anhf)Y lm {f) 



exp(-p)P JV (p)y LM (p) (Al) 



Note f = r/r max and p = p/T p . If r max and/or T p vary 
with time, then one needs to include additional terms 
which we do not elaborate in details. To obtain the evo- 
lution equation for the expansion coefficients a^^ M , we 
define the following shorthand to project out the expan- 
sion coefficients from any function F, 

(n'l'm'N'L'M',F(x,p)) = / f 2 dr dfl f dp- 
Jo J Jo 

dn p R nn , (<*„»,» , f)y,r m , (?)p N , (p)Y£, M , (p)f(*, p ) (A2) 



Performing such projection for the first term, we obtain 
(nWiVW,M|M } = ^p (A3 ) 



The use of spherical basis vectors is convenient as three 
components of a vector V arc directly related to spherical 
harmonics Yy, 



An 

v t = |v|W yF M (y) 



(A6) 



The gradient formula Eq. (|A4j) can be further simplified 
if one has spherical Bcsscl function, F(r) = ji(kr), with 
the help of the following recurrence relations 

—jl(kr) = kji-x(kr) - l —^-ji(kr) 
dr r 

^ji(kr) = -kj l+1 (kr) + -ji(kr) (A7) 
dr r 



Applying to the first and second terms in Eq. (|A4[) . one 
has for our case 



Vji(a n if)Y lm (r) 



y 2(21 + iy J V( a "i f ) y r,m-i(r)G, m - i, l,i\l,m) (A8) 

where we have combined two terms together into a com- 
pact form. The second term becomes 



NLMt.\ a nl 



W/(x,P,t)= J2 a nL M {t) 



V2 



nlmNLM 



^max 



E X^M+i + 5 T,l-l)J + ^ M^nir)Y Im ^(r) 



i 1 



2(2/ + 1)- 



/4tt, 



dt 



dt 



(l,m-i, l,i\l,m)exp(-p)P N (p)Y LM (p)d yYii(p) (A9) 

Following the same procedure as done for the first term, 
we project out the expansion coefficients a^h^ f° r the 
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second term, 

(n'l'm'N'L'M',vVm,p,t))= £ ^>™(i) 



nlmNLM i 



aw 



z + z + i 



^ 2(2? + 1) 
V2 ^ 



(Z, m — i, 1, to) 



f drj v {a n , v r)ji(a n ir) 



ji'+i(a n 'i') ji+i{a n i) Jo 

poo 

!Oyj? m ,(f)y r , m _,(r) y j5 2 dpP^(p)exp(-p)/V(p) 

/47T 

dn p YZ, M ,(p)Y LM (p)\j y y ^(p) ( A1 °) 



The integrals J e?£7 and J dp can be done using orthogonal 
relations for spherical harmonics and Laguerre function. 
The integral J d£l p involves the product of three spherical 
harmonics, which can performed with the help of the 
following relation, 



dn p Yl, MI {p)Y LM {p)Yiiij>) 



3(2i+l) 



/ 47r(2i' + 1) 
(L, 0, 1, 0\L', 0)(L, M, 1, i\L', M') (All) 



Note that the above Clebsch-Gordan coefficients are non- 
zero only when \L — L'\ < 1 and L — V — 1 are even 
numbers. This implies that L = V ± 1 and M = M' — i. 
The final result for the second term is, 



i'l'm'N>L>M',v.Vf(Z,p,t))=J2 E a ™(*) 



nlmNLM i 



a n i 



I r [n! , I', n, 1}6 nn >(5l>j +1 + dVj-i)i 



ll + V + 1 
2(2/ + 1) 



(6l>,l+i + 6 L ',l-i)J^j^(L, 0, 1, 0\L', 0) 
S m ', m -iQ', m - i, 1, i\l, m)6 M ',M+i(L, M, 1, i\L', M') 



(A12) 



where I r [n' , V , n, I] has been defined in Eq. (j2T))) . Combin- 
ing with the first term, we finish the multipole expansion 
of the left hand side in Boltzmann equation. 
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